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Q This paper deals with stabihty of a certain class of fractional order linear and nonlinear 

l/^ systems. The stability is investigated in the time domain and the frequency domain. The 

general stability conditions and several illustrative examples are presented as well. 
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I 

g 1 Introduction 

^ Fractional calculus is more than 300 years old topic. A number of applications where fractional 
^ calculus has been used rapidly grows. This mathematical phenomena allow to describe a real 
O object more accurate that the classical "integer" methods. The real objects arc generally 
^ fractional [40,42,50,70,72], however, for many of them the fractionality is very low. The main 

^ J reason for using the integer-order models was the absence of solution methods for fractional 

^ differential equations. 

Q Recently, the fractional order hnear time invariant (FOLTI) systems have attracted lots 

^ of attention in control systems society (e.g.: [16,34,42,52,54]) even though fractional-order 

^ control problems were investigated as early as 1960s [35]. In the fractional order controller, 
^ the fractional order integration or derivative of the output error is used for the current control 
force calculation. 

The fractional order calculus plays an important role in physics [45,60,67], thermodynamics 
[31,57], electrical circuits theory and fractances [7,12,14,21,40,71], mechatronics systems [55], 
signal processing [56,68], chemical mixing [41], chaos theory [62,64], and biological system 
as well [23]. It is recommended to refer to (e.g.: [6,39,43,59,73]) for the further engineering 
applications of fractional order systems. The question of stability is very important especially in 
control theory. In the field of fractional-order control systems, there are many challenging and 
unsolved problems related to stability theory such as robust stability, bounded input - bounded 
output stability, internal stability, root-locus, robust controllability, robust observability, etc. 
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For distributed parameter systems with a distributed delay [44], provided an stability anal- 
ysis method which may be used to test the stability of fractional order differential equations. 
In [13], the co-prime factorization method is used for stability analysis of fractional differential 
systems. In [36] , the stability conditions for commensurate FOLTI system have been provided. 
However, the general robust stability test procedure and proof of the validity for the general 
type of the FOLTI system is still open and discussed in [48] . Stabihty has also been investigated 
for fractional order nonlinear system (chaotic system) with commensurate and incomensurate 
order as well [2,62,63]. 

This paper is organized as follows. In Sec. 2 is briefly introduced the fractional calculus. 
Sec. 3 is on fractional order systems. In Sec. 4 are analyzed the stability conditions of fractional 
order linear and nonlinear systems. Sec. 5 concludes this paper with some remarks. 



2 Fractional Calculus Fundamentals 

2.1 Definitions of Fractional Derivatives and Integrals 

The idea of fractional calculus has been known since the development of the regular calculus, 
with the flrst reference probably being associated with Leibniz and L'Hospital in 1695 where 
half-order derivative was mentioned. 

Fractional calculus is a generalization of integration and differentiation to non-integer order 
fundamental operator a-D[, where a and t are the limits of the operation and r & R. The 
continuous integro-differential operator is defined as 

f : r > 0, 

aDl={ 1 :r = 0, 

fjdr)-^ :r<0. 

The three definitions used for the general fractional differintegral are the Grunwald-Letnikov 
(GL) definition, the Riemann-Liouville (RL) and the Caputo definition [41,50]. The GL is given 
here 



[— ] 

aDlfit) = lim h-^ Yl (-1)' ( •) /(^ - ^■^)' 

j=Q V-?/ 



(1) 



where [.] means the integer part. The RL definition is given as 



for (n — 1 < r < n) and where r(.) is the Gamma function. The Caputo's definition can be 
written as 



1 r f 



n). 



T 



for (n — 1 < r < n). The initial conditions for the fractional order differential equations with 
the Caputo's derivatives are in the same form as for the integer-order differential equations. 



2 



2.2 Some Properties of Fractional Derivatives and Integrals 

The main properties of fractional derivatives and integrals are the following: 

1. If f{t) is an analytical function of t, then its fractional derivative QD"f{t) is an analytical 
function of t, a. 

2. For a = n, where n is integer, the operation qD^ f{t) gives the same result as classical 
differentiation of integer order n. 

3. For a = the operation QDff{t) is the identity operator: 

oAVW = fit) 

4. Fractional differentiation and fractional integration are linear operations: 

aDl{\f{t)+^^g{t)) = XaDlf{t)+f,,Dlg{t). 

5. The additive index law (semigroup property) 

oD%D?f{t) = oD%D^f{t) = oDr^fit) 

holds under some reasonable constraints on the function f{t). 

The fractional-order derivative commutes with integer-order derivation 



under the condition t = a we have f^''\a) = 0, (A; = 0, 1, 2, . . . , n — 1). The relationship 
above says the operators ^ and commute. 

6. The formula for the Laplace transform of the RL fractional derivative ^ has the form [50] : 

/•oo n— 1 

/ e-^* oDUit) dt = s^F{s) A^-'- V(t) L=o ' 



k=0 

for (n— 1 < r < n), where s = juj denotes the Laplace operator. For zero initial conditions, 
Laplace transform of fractional derivatives (Grunwald-Letnikov, Riemann-Liouville, and 
Caputo's), reduces to: 

L{,Dlf{t)} = s^F{s). 

7. Geometric and physical interpretation of fractional integration and fractional differentia- 
tion were exactly described in Podlubny's work [51]. 

Some others important properties of the fractional derivatives and integrals as for example 
Leibniz's rule, translation. Chain rule, bahaviour and dependence on limit and so on, we can 
find out in several works (e.g.: [41,42,50], etc.). 
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3 Fractional- Order Systems 



3.1 Fractional LTI Systems 

A general fractional-order system can be described by a fractional differential equation of the 
form 



= h,nD^"'u{t) + hra-iD^'^-^u{t) + . . . + boD^"u{t), 



(4) 



or by the corresponding transfer function of incommensurate real orders of the following form 
[50]: 

^ bms^"^ + ... + bis^^ + fepg^" ^ Q{s^^) 
anS""" + ... + ais°i + aos"o P(s"fe) ' 
where D'^ = qD] denotes the Riemann-Liouville or Caputo fractional derivative [50]; {k = 
0, . . . n), bk {k = 0, . . . m) are constant; and {k = 0, . . . n), Pk {k = 0, . . . m) are arbitrary 
real numbers and without loss of generality they can be arranged as q;„ > q;„_i > . . . > ao, 
and Prn > Pm-i > ...> Po- 

The incommensurate order system ^ can also be expressed in commensurate form by the 
multi-valued transfer function [9] 

+ ■ ■ ■ + feis^/^ + bo 



His) 



(v > 1). 



(6) 



+ aisV^ + ao 

Note that every fractional order system can be expressed in the form ^ and domain of the 
H[s) definition is a Riemann surface with v Riemann sheets [32]. 

In the particular case of commensurate order systems, it holds that, = ak, j3k = ak, (0 < 
a < l),Vfc G Z, and the transfer function has the following form: 



G{s) = 



E 



(7) 



.k=0 (^k[s^ 

With N > M, the function G{s) becomes a proper rational function in the complex variable 
which can be expanded in partial fractions of the following form: 



G{s) = Ko 



■ N 

E 

i=l 



A, 



where Aj [i = 1, 2, .., A^) are the roots of the pseudo-polynomial P(s°) or the system poles which 
are assumed to be simple without loss of generality. The analytical solution of the system ([8 
can be expressed as 



yit) = { 



■ N 

i=l 



A, 



s" + Ai 



(9) 



A fractional order plant to be controlled can be described by a typical ra-term linear homo- 
geneous fractional order differential equation (FODE) in time domain 



a„ D'i-y{t) + ■ ■ ■ + ai D^yit) + ao DTyit) = 



(10) 
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where ak{k = 0, 1, ■ ■ ■ , n) are constant coefficients of the FODE; a^, (fc = 0, 1, 2, ■ ■ ■ , n) are real 
numbers. Without loss of generality, assume that a„ > > ■ ■ ■ > cto > 0. 
The analytical solution of the FODE (10) is given by general formula [50] 

X °° (—1)"* 
an ^-^ ml ^-^ 

m=0 kQ + ki + ... + k„_2=m 

feO>0;... ,fc„_2>0 

n-2 / ^ \ ki 

X 



i=0 
n-2 

+ + (11) 

j=0 

where (m; fco, ki, . . . , A;„_2) are the multinomial coefficients and Sk{t,y; fi, v) is the function of 
Mittag-Leffier type introduced by Podlubny [50]. The function is defined by 

£,{t,y-ii,u)=t^^^^-'Efl{yt^), (A; = 0, 1, 2, . . .), (12) 

where E^^y{z) is the Mittag-Leffier function of two parameters [27]: 

00 j 



i=0 



where e.g. Ei^i{z) = e^, and where its fc-th derivative is given by 



Consider a control function which acts on the FODE system (10) as follows: 

an Dt-y{t) + ■ ■ ■ + ai D^^y{t) + a^ Dt"y{t) = u{t). (15) 
By Laplace transform, we can get a fractional transfer function: 

G{s) = ^ = ^ . (16) 

The fractional order linear time-invariant system can also be represented by the following state- 
space model 

oD^x{t) = Ax{t) + Bu{t) 

y{t) = Cx{t) (17) 

where x G -R", u G K'' and y & BP are the state, input and output vectors of the system 
and A G i?"^", B G i?"'^'', C G i?^^", and q = [gi, g2, • • • , (Jn]^ are the fractional orders. 
If qi = q2 = ■■■Qn, system (17) is called a commensurate order system, otherwise it is an 
incommensurate order system. 
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A fractional-order system described by n-term fractional differential equation (15) can be 
rewritten to the state-space representation in the form [17,69]: 










1 . . " 




Xi(t) 














1 . 




X2{t) 



















+ 










— ai/a„ . . a„_i/a„ 




Xn{t) 







u{t) 



[I 



] 



Xl{t) 
X2{t) 



Xn{t) 



(18) 



where ao = 0, gi = ai, g2 = «n-i — «n-2, ■ ■ ■ In = Oin — On-i, and with initial conditions: 

xi(0) = XQ^=yQ, X2(0) = Xq^^ = 0, . . . 

» _ J yf\ if « = 2/i; + 1, 



x.(0) 



X 



0, if i = 2A;, 



i < n. 



(19) 



The n-term FODE (15) is equivalent to the system of equations (18) with the initial condi 



tions (19). 



Similar to conventional observability and controllability concept, the controllability is de- 



fined as follow [38]: System (17) is controllable on [to^/maz] if controllability matrix Ca 



[BlABlA^Bl . . . 1^4" ^B] has rank n. The observability is defined as follow [38]: System (17) is 
observable on [to,^/mai] if observanility matrix Oa = [C|CA|CA^| . . . |Cy4"'~-^]'^ has rank n. 

3.2 Fractional Nonlinear Systems 

Generally, we consider the following incommensurate fractional order nonlinear system in the 
form: 



0DfXi{t) = fi{Xi{t),X2(t),...,Xn{t),t) 



Xi{0) 



where q are initial conditins, or in its vector representation: 

D'lx = f(x), 

where q = [gi, q2, ■ ■ ■ , qn]"^ for < < 2, (i = 1, 2, . . . , n) and x G -R". 



(20) 



(21) 



The equilibrium points of system (21) are calculated via solving the following equation 

f(x) = (22) 



and we suppose that equilibrium point of system (plj). 
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4 Stability of the Fractional Order Systems 



4.1 Preliminary Consideration 

Stability as an extremely important property of the dynamical systems can be investigated in 
various domain [18,19]. Usual concept of bounded input - bounded output (BIBO) or external 
stability in time domain can be defined via the following general stability conditions [37]: 

A causal LTI system with impulse response h{t) to be BIBO stable if the necessary and 
sufficient condition is satisfied 

/■oo 

\\h{T)\\dT < oo. 







where output of the system is defined by convolution 

POO 

y[t) = h{t) * U{t) = / h{T)u{t - T)dT, 



where u,y E L^o and h E Li. 

Another very important domain is frequency domain. In the case of frequency method 
for evaluating the stability we transform the s-plane into the complex plane Go{juj) and the 
transformation is realized according to the transfer function of the open loop system Go{juj). 
During the transformation, all roots of the characteristic polynomial are mapped from s-plane 
into the critical point (— 1, jO) in the plane Go{juj). The mapping of the s-plane into Go{juj) 
plane is conformal, that is, the direction and location of points in the s-plane is preserved in 
the Go{juj) plane. Frequency investigation method and utilization of the Nyquist frequency 
characteristics based on argument principle were described in the paper [46]. 

However, we can not directly use an algebraic tools as for example Routh-Hurwitz criteria 
for the fractional order system because we do not have a characteristic polynomial but pseudo- 
polynomial with rational power - multivalued function. It is possible only in some special 
cases [2]. Moreover, modern control method as for example LMI (Linear Matrix Inequality) 
methods [43] or other algorithms [29, 30] already have been developed. The advantage of 
LMI methods in control theory is due their connection with the Lyapunov method (existence 
a quadratic Lyapunov function). More generally, LMI methods are useful to test of matrix 
eigenvalues belong to a certain region in complex plane. A simple test can be used [3]. Roots of 
polynomial -P(s) = det (sJ — A) lie inside in region —it/ 2 — 6 < arg(s) < 7t/2 + S if eigenvalues 
of the matrix 



A, 



A cos 6 —A sin 6 
A sin 6 A cos 6 



= A 



cos 6 — sin (5 
sin S cos 6 



(23) 



have negative real part, where ® denotes Kronecker product. This property has been used to 
stability analysis of ordinary fractional order LTI system and also for interval fractional order 
LTI system [65]. 

When dealing with incommensurate fractional order systems (or, in general, with fractional 
order systems) it is important to bear in mind that P(s"), a G R is a multivalued function 
of s°, a = -, the domain of which can be viewed as a Riemann surface with finite number of 
Riemann sheets v, where origin is a branch point and the branch cut is assumed at R~ (see 
Fig. [T]). Function s° becomes holomorphic in the complement of the branch cut line. It is 
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a fact that in multivalued functions only the first Riemann sheet has its physical significance 
[28]. Note that each Riemann sheet has only one edge at branch cut and not only poles and 
singularities originated from the characteristic equation, but branch points and branch cut of 
given multivalued functions are also important for the stability analysis [10]. 



Im 




Figure 1: Branch cut (0, — oo) for branch points in the complex plane. 

In this paper the branch cut is assumed at R^ and the first Riemann sheet is denoted by Q 
and defined as 

Q := {re^^\r > 0,-7i< (p<n}. (24) 

It is well-known that an integer order LTI system is stable if all the roots of the characteristic 
polynomial P{s) are negative or have negative real parts if they are complex conjugate (e.g.: 
[18]). This means that they are located on the left of the imaginary axis of the complex s-plane. 
System G{s) = Q{s)/P{s) is BIBO stable if 

3, 11^(5)11 < M < oo, M>0, Vs,3fJ(s)>0. 

A necessary and sufficient condition for the asymptotic stability is [25]: 

limt^^\\X{t)\\=0. 

According the final value theorem proposed in [26], for fractional order case, when there is 
a branch point at s = 0, we assume that G{s) is multivalued function of s, then 

x{oo) = lims^Q[sG{s)]. 

Example 1: Let us investigate the simplest multi- valued function defined as follow 

w = (25) 

and there will be two s-planes which map onto a single w-plane. The interpretation of the two 
sheets of the Riemann surface and the branch cut is depicted in Fig. [2j 
Define the principal square root function as 

111 

fi{s) = |s| 2e 2 = re 2 , 
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Figure 2: Riemann surface interpretation of the function w = . 

where r > and — tt < < +7r. The function is a branch of w. Using the same notation, 
we can find other branches of the square root function. For example, if we let 

J2[s) = |s| 2e 2 = re 2 , 

then /2(s) = — and it can be thought of as "plus" and "minus" square root functions. The 
negative real axis is called a branch cut for the functions fi{s) and /2(s). Each point on the 
branch cut is a point of discontinuity for both functions fi{s) and /2(s). As has been shown 



in [32], the function described by (25) has a branch point of order 1 at s = and at infinity. 
They are located at ends of the branch cut (see also Fig. [T|. 

Example 2: Let us investigate the transfer function of fractional-order system (multivalued 
function) defined as 

Gis) = (26) 
where a e R {0 < a < 2) and b e R {b> 0). 



The analytical solution of the fractional order system (26) obtained according to relation 



(11) has the following form: 

g{t) = So{t,-b;a,a). (27) 



The Riemann surface of the function ( 26 ) contains an infinite number of sheets and infinitely 
many poles in positions 

1 j(7r + 27rn) 

s = b^e — 5 — , n = 0,±1,±2,. . . , for (a > 0) and (6 > 0). 

The sheets of the Riemann surface are all different if a is irrational. 

For 1 < a < 2 we have two poles corresponding to n = and n = —1, and poles are 

s = b^e a . 



However, for < a < 1 in (26) the denominator is a multivalued function and singularity 



of system can not be defined unless it is made singlevalued. Therefore we will use the Riemann 
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surface. Let us investigate transfer function (26) for a = 0.5 (half-order system), then we get 



Gis) = (28) 

and by equating the denominator to zero we have 

si + 6 = 0. 

Rewriting the complex operator in exponential form and using the well known relation 
e^'^ + 1 = (or _)_ i = g) we get the following formula: 



From relationship (29) can be deduced that the modulus and phase (arg) of the pole are: 

r = b^ and = ±2n{l + k) for A; = 0, 1, 2, . . . 

However the first sheet of the Riemann surface is defined for range of — tt < (j) < +7i, the pole 
with the angle = ±2tt does not fall within this range but pole with the angle = 27r falls 
to the range of the second sheet defined for vr < < 27r. Therefore this half-order pole with 
magnitude 6^ is located on the second sheet of the Riemann surface that consequently maps 
to the left side of the w-plane (see Fig. [s]). On this plane the magnitude and phase of the 
singlevalued pole are 6^ and vr, respectively [32]. 

Example 3: Analogous to previous examples we can also investigate function 

U7 = S5, (30) 

where in this case the Riemann surface has three sheets and each maps onto one-third of the 
w-plane (see Fig. [i]). 

Definition 1. Generally, for the multivalued function defined as follow 

w = s"" , (31) 

where v G (f = 1, 2, 3, . . . ) we get the v sheets in the Riemann surface. In Figjs] is shown 
the relationship between the w-plane and the v sheets of the Riemann surface where sector 
—ir/v < arg(w) < n/v corresponds to Q (first Riemann sheet). 

Definition 2. Mapping the poles from s^-plane into the w-plane, where q & Q such as g = ^ 
for k,m & N and |arg(w)| = |0|, can be done by the following rule: If we assume k = 1, then the 
mapping from s-plane to w-plane is independent of k. Unstable region from s-plane transforms 
to sector \4>\ < ^ and stable region transforms to sector ^ < |0| < ^. The region where 
101 > ^ is not physical. Therefore, the system will be stable if all roots in the ly-plane lie in the 
region |0| > Stability regions depicted in Fig. [6] correspond to the following propositions: 



1. For k < m {q < 1) the stability region is depicted in Fig. 6 (a 



2. For k = m (g = 1) the stability region corresponds to the s-plane (see Fig. |6(b) ). 



3. For k > m (g > 1) the stability region is depicted in Fig. 6(c 
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Figure 3: Correspondence between the s-plane and the w-plane for Eq.(28) 



4.2 Stability of Fractional LTI Systems 

As we can see in previous subsection, in the fractional case, the stabihty is different from the 
integer one. Interesting notion is that a stable fractional system may have roots in right half 
of complex w-plane (see Fig. [6]). 

Since the principal sheet of the Riemann surface is defined —it < arg(s) < tt, by using the 
mapping w = s'^, the corresponding w domain is defined by —qn < a,Yg{w) < qir, and the w plane 
region corresponding to the right half plane of this sheet is defined by — g7r/2 < arg(w) < g7r/2. 

Consider the fractional order pseudo-polynomial 

Q{s) = ais'^' + a2S^2 + . . . + a^s"?" = ais^'/'^' + ass"^/'^^ + . . . + a„s'="/'^", 

where are rational number expressed as Ci/di and are the real numbers for i = 1,2, . . . ,n. 
If for some i, Ci = then di = 1. Let v be the least common multiple (LCM) of di,d2, . . . dn 
denote as f = LCM{(ii, d2, . . . d^}, then [26] 

Q(s) = ais^ + a2S^ + . . . + a^s"^ = aiis'-Y' + 02(5^)^^ + . . . + a^is-Y". (32) 

The fractional degree (FDEG) of the polynomial Q{s) is defined as [26] 

FDEG{(5(s)} = max{'Ui, t;2, . . . , w„}. 
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Im(w) 




Range for 
sheet 1 



Re(w) 



(a) Riemann surface 



(b) Complex w-plane 



Figure 4: Correspondence between the 3-sheets Riemann surface and w-plane for Eq.(30) 



Im(w) 




sheet 2 




> Re(w) 



sheet V 



(a) Riemann surface 



(b) Complex w-plane 



Figure 5: Correspondence between the w-plane and the Riemann sheets for Eq.(31). 



The domain of definition for (32) is the Riemann surface with v Riemann sheets where origin 



is a branch point of order f — 1 and the branch cut is assumed at R . Number of roots for 



fractional algebraic equation (32) is given by the following proposition [8]: 



Proposition 1. Let Q{s) be a fractional order polynomial with FDEG{Q(s)} = n. Then the 
equation Q(s)=0 has exactly n roots on the Riemann surface [8]. 
Definition 3. The fractional order polynomial 



Q[s) = aiS" + « + . . . + a„S" + a„+i 
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Irrii 



stable 
region 



n 



Re 

unstable stable 
region region 




(a) < q < 1 



(b) <? = 1 



(c) 1 < 9 < 2 



Figure 6: Stability regions of the fractional order system. 



is minimal if FDEG{(5(s)} = n. We will assume that all fractional order polynomial are 
minimal. This ensures that there is no redundancy in the number of the Riemann sheets [26]. 

On the other hand, it has been shown, by several authors and by using several methods, that 
for the case of FOLTI system of commensurate order, a geometrical method of complex analysis 
based on the argument principle of the roots of the characteristic equation (a polynomial in 
this particular case) can be used for the stability check in the BIBO sense (see e.g. [37,46]). 
The stability condition can then be stated as follows [36,37,58]: 

Theorem 1. A commensurate order system described by a rational transfer function ([T]) is 
stable if only if 

vr 

|arg(Aj)| > a — , for all i 

with Aj the i-th root of Pi^s"^). 

For the FOLTI system with commensurate order where the system poles are in general 
complex conjugate, the stability condition can also be expressed as follows [36,37]: 
Theorem 2. A commensurate order system described by a rational transfer function 



G{w) 



OH 

PivS] 



where w = s"^, q E R"*", (0 < g < 2), is stable if only if 

|arg(wi)| > g^, 

with Wwi G C the i-th root of P{w) = 0. 

When w = is a single root (singularity at the origin) of P, the system cannot be stable. 
For g = 1, this is the classical theorem of pole location in the complex plane: have no pole 
in the closed right half plane of the first Riemann sheet. The stability region suggested by 
this theorem tends to the whole s-plane when q tends to 0, corresponds to the Routh-Hurwitz 
stability when q = 1, and tends to the negative real axis when q tends to 2. 
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Theorem 3. It has been shown that commensurate system (17) is stable if the following 

(33) 



condition is satisfied (also if the triplet A, B, C is minimal) [4,37,61-63]: 



|arg(eig(A))| > q 



IT 



where < g < 2 and eig(A) represents the eigenvalues of matrix A. 

Proposition 2. We can assume, that some incommensurate order systems described by the 



FODE (15) or (17), can be decomposed to the following modal form of the fractional transfer 



function (so called Laguerre functions [5]): 

F{s] 



A, 



(34) 



i=l k=l 

for some complex numbers Aj^^, Aj, and positive integer rzfc. 



A system (34) is BIBO stable if and only if qi and the argument of A, denoted by arg(Aj 



in (34) satisfy the inequalities 

< gj < 2 and |arg (Aj)| < vr ( 1 



for all i. 



(35) 



Henceforth, we will restrict the parameters to the interval G (0, 2). For the case g^ = 1 for 
all i we obtain a classical stability condition for integer order system (no pole is in right half 



plane). The inequalities (35) were obtained by applying the stability results given in [1,37]. 
Theorem 4. Consider the following autonomous system for internal stability definition [15]: 



Xo, 



(36) 



with q = [gi, g2, . . . , qn]"^ and its n-dimensional representation: 



oD't'xiit) 
oDfx^it) 



aiiXi(t) + ai2X2{t) H h ainXnit) 

a2lXi{t) + a22X2{t) H h a2nXn{t) 



oDl"Xn(t) = a„iXi(t) + an2X2{t) H h annXn{t) 



(37) 



where all g^'s are rational numbers between and 2. Assume m be the LCM of the denominators 
Tij's of gj's , where qi = Vi/ui, Vi, Ui G for i = 1, 2, . . . , n and we set 7 = 1/m. Define: 



/A 



det 



mqi 



— an —<Xi2 
a2i A"^92 - a22 



\ 



— Q-la 
—a-2n 

A*""" - ar, 



\ 



J 



(38) 



The characteristic equation (38) can be transformed to integer order polynomial equation if 



all gi's are rational number. Then the zero solution of system (37) is globally asymptotically 



stable if all roots Aj's of the characteristic (polynomial) equation (38) satisfy 



TT 



|arg(Ai)| > 7- for all i. 
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Denote A by s'^ in equation (38 ), we get the characteristic equation in the form det{s'^ I — A) = 
and this assumption was proved in paper [15]. 

Corollary 1. Suppose qi = q2 = ■ ■ ■ ,(ln = Q, Q ^ (0,2), all eigenvalues A of matrix A in 
(18) satisfy |arg(A)| > q7i/2, the characteristic equation becomes det(s'^J — A) = and all 
characteristic roots of the system (17) have negative real parts [15]. This result is Theorem 1 
of paper [36]. 

Remark 1. Generally, when we assume s = Irle**^, where |r| is modulus and (p is argument of 
complex number in s-plane, respectively, transformation w = to complex w-plane can be 
viewed as s = |r|mem and thus |arg(s)| = m.|arg(w)| and \s\ = Proof of this statement 

is obvious. 

Stability analysis criteria for a general FOLTI system can be summarized as follow: 

The characteristic equation of a general LTI fractional order system of the form: 



a„s"" + . . . + ais"i + aos"° = ^ a^s"' = (39) 

4=0 

may be rewritten as 



n 
i=0 



and transformed into w-plane 



ttiw' = 0, (40) 



i=0 



with w = S"i^ where m is the LCM of fj. The procedure of stability analysis is (see e.g. [53]): 



1. For given Oj calculate the roots of Eq.(40) and find the absolute phase of all roots 

2. Roots in the primary sheet of the w-plane which have corresponding roots in the s-plane 
can be obtained by finding all roots which lie in the region \(f)w\ < ^ then applying the 
inverse transformation s = w"^ (see Remark 1.). The region where 10^^,1 > ^ is not 
physical. For testing the roots in desired region the matrix approach can be used (23). 

3. The condition for stability is ^ < |0^| < ^. Condition for oscillation is |0^| = ^ 
otherwise the system is unstable (see Fig. 5(b)). If there is not root in the physical 
s-plane, the system will always be stable [53]. 

Example 4. Let us consider the linear fractional order LTI system described by the transfer 
function [16,50]: 

Gis) = ^ = 1 f41^ 

^' U{s) 0.8s2-2 + 0.5s0-9 + l' ^ ^ 

and corresponding FODE has the following form: 

0.8 oDry{t) + 0.5 ^D^^yit) + y{t) = u{t) (42) 
with zero initial conditions. 
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The system (42) can be rewritten to its state space representation (xi(t) = y(t)): 



oDioxi{t) 



y{t) 



1 
-1/0.8 -0.5/0 

Xi{t) 



xi{t) 

X2{t) 



+ 





1/0.^ 



[1 0] 



X2{t) 



The eigenvalues of the matrix A are Ai,2 = —0.3125 ± 1.0735j and then |arg(Ai,2)| 



(43) 
1.8541. 



Because of various derivative orders in (43), the Theorem 3 cannot be used directly. 




Figure 7: Analytical solution of the FODE (42) where u{t) = for 50 sec. 



The analytical solution of the FODE (42) for u{t) = obtained from general solution (11) 
has form: 



1 >^ (-i)^ /2. 

08 ^ k\ V 

fc=0 ^ 



£,{t,-^;2.2- 0.9,2.2 + 0.9k). 
U.o 



(44) 



In Fig. u^is depicted the analytical solution of the FODE (42) where u{t) = 0. As we can see 



in the figure, solution is stable because limt^^y{t) = 0. Let us investigate stability according 
to the previously described method. The corresponding characteristic equation of system is: 



P(s) : 0.8s2-2 + 0.5s°-^ + l = 0.8s w + 0.5s 1 = 0, (45) 
when m = 10, w = sra then the roots Wj's and their appropriate arguments of polynomial 



P{w) : 0.8w^^ + 0.5w 







(46) 



are: 



Wl,2 



-0.9970 ±0.1182j, |arg(w;i,2)| = 3.023; ^3,4 = -0.9297 ± 0.4414j, |arg(«;3,4)| = 2.698; 
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W5,6 = -0.7465 ± 0.6420J, |arg(w;5,6)| 
W9,io = -0.259 ±0.9625j, |arg(w;9,io)| 
wi^^ii = 0.3080 ± 0.9772j, |arg(wii,i2) 
Wi7,i8 = 0.7793 ± 0.6795J, |arg(w;i7,i8) 
W2i,22 = 1.0045 ± 0.1684J, |arg(w2i,22) 



2.431; w;7,8 = ■ 
1.834; w;n,i2 = 
= 1.265; Wi5,i6 
= 0.717; Wi9,2o 
= 0.1661; 



-0.5661 ± 0.8633J, |arg(w7,8)| = 2.151; 
-0.0254±1.0111j, |arg(wii,i2)| = 1.595 
= 0.5243 ±0.8359j, |arg(«;i5,i6)| = 1.010 
= 0.9084 ±0.3960j, |arg(«;i9,2o)| = 0.411 



Physical significance roots are in the first Riemann sheet, which is expressed by relation 
—Tc/m < (f) < tt/itl, where (f) = aig{w). In this case they are complex conjugate roots tL'2i,22 = 
1.0045±0.1684j (|arg(w2i,22)| = 0.1661), which satisfy conditions |arg(i(72i,22)| > vr/2m = 7r/20. 



It means that system (42) is stable (see Fig. |8]). Other roots of the polynomial equation (46) 
lie in region |0| > — which is not physical (outside of closed angular sector limited by thick 



line in Fig. 8(b) ). 





(a) 10-sheets Riemann surface 



(b) Poles in complex w-plane 



Figure 8: Riemann surface of function w = sw and roots of Eq.(46) in complex ly-plane 



In Fig. 



is depicted the Riemann surface of the function w = sw with the 10-Riemann 



sheets and in Fig. 8(b) are depicted the roots in complex w-plane with angular sector corre- 



sponds to stability region (dashed line) and the first Riemann sheet (thick line). 

The interesting notion of Remark 1 should be mentioned here. The characteristic equation 



(45) has the following poles: 



si 2 = -0.10841 ± 1.19699j, 

in the first Riemann sheet in s-plane, which can be obtained e.g. via the Matlab routine as for 
instance: 

»s=solve('0.8*s"2. 2+0. 5*s~0. 9+1=0' , 's') 

When we compare |arg(t(;2i,22)| = 0.1661 and |arg(si^2)| = 1-661, we can see that |arg(si_2)| = 
m|arg(w2i,22)|; where m = 10 in transformation w = . The first Riemann sheet is trans- 
formed from s-plane to w-plane as follow: — tt/IO < arg(w) < vr/lO and in order to —it < 
lO.arg(w) < vr. Therefore from this consideration we then obtain |arg(s)| = 10.|arg(w)|. 
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Example 5. Let us examine an interesting example of application, so called Bessel function 
of the first kind, which transfer function is [37]: 

His) = v., ^{s) > 0. (47) 

Vs^ + 1 

We have two branch points si = i, and S2 = —i and two cuts. One along the half line {—oo+i, i) 
and another one along the half line (— oo — i, —i). In this doubly cut complex plane, we have 
the identity + 1 = ^/s — iy/s + i. The well known asymptotic expansion of Eq.(47) is: 

Mt)-/|cos(t-^) = yfHE,, (-(t-J 

According to the branch points and above asymptotic expansion we can state, that system 
described by the Bessel function (47) is on boundary of stability and has oscillation behaviour. 
Example 6. Consider the closed loop system with controlled system (electrical heater) 

^<^' = 39.96..^ + 0.598 ^''^ 

and fractional order controller 

C(s) = 64.47+ 12.46s (49) 
The resulting closed loop transfer function Gc{s) becomes [49]: 

G(s) = ^= 12.46. + 64.47 

^ W{s) 39.69si-25 + 12.46s + 65.068 ^ ' 



The analytical solution (impulse response) of the fractional order control system ( 50 ) is: 

\k / 1 O 



12.46 -ir /12.46\ 65.068 

y(t) = x^Jt, : 1.25, 0.25 -A; 

^ 39.69 ^ k\ 1 39.69/ ^' 39.69' ' ^ 

64.47 ^(-1)*^ /65.068\'' ^, 12.46 

+ > x^fe t, ; 1.25 - 1,1.25 + A; 51 

39.69 ^ k\ V 39.69 / ^' 39.69' ' ' ^ ' 



fc=0 

with zero initial conditions. 

The characteristic equation of this system is 

39.69s^-2^ + 12.46s + 65.068 = ^ 39.69st + 12.46st + 65.068 = (52) 
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Using the notation = s™, where LCM is m = 4, we obtain a polynomial of complex variable 
w in form 

39.69m;^ + 12.46ti;^ + 65.068 = 0. (53) 



Solving the polynomial (53) we get the following roots and their arguments: 

Wx = —1.17474, |arg(wi)| = n 

W2,3 = -0.40540 ± 1.0426J, |arg(w2,3)| = 1.9416 

m;4,5 = 0.83580 ± 0.64536j, |arg(^^;4,5)| = 0.6575 

This first Riemann sheet is defined as a sector in ly-plane within interval —tt/4 < aig{w) < tt/4. 
Complex conjugate roots ^4^5 lie in this interval and satisfies the stability condition given as 
|arg(w)| > |, therefore system is stable. The region where |arg(w)| > | is not physical. 
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4.3 Stability of Fractional Nonlinear Systems 

As it was mentioned in [36], exponential stability cannot be used to characterize asymptotic 
stability of fractional order systems. A new definition was introduced [43]. 



Definition 3. Trajectory x{t) = of the system (20) is t ^ asymptotically stable if there is 
a positive real q such that: 

V||x(t)|| with t< to, 3 N{x{t)), such that > to,\\x{t)\\ < Nt~'^. 

The fact that the components of x{t) slowly decay towards following leads to fractional 
systems sometimes being called long memory systems. Power law stability t~'^ is a special case 
of the Mittag-Leffler stability [33]. 

According to stability theorem defined in [66], the equilibrium points are asymptotically 
stable for gi = g2 = • • • = Q'n = <? if all the eigenvalues Aj, {i = 1,2, ... ,n) of the Jacobian 
matrix J = 9f/9x, where f = [/i, /2, /n]^, evaluated at the equilibrium, satisfy the 

condition [61,62]: 

|arg(eig(J))| = |arg(Ai)| > g-, 2 = 1,2, ...,n. (54) 

Fig. [6] shows stable and unstable regions of the complex plane for such case. 

Now, consider the incommensurate fractional order system qi 7^ q2 Qn and suppose 

that m is the LCM of the denominators -Uj's of g^'s, where g^ = Vi/ui, Vi,Ui E for i = 



1,2, ... ,n and we set 7 = 1/m. System (21 ) is asymptotically stable if: 

|arg(A)| > 7^ 

for all roots A of the following equation 

det(diag([A'"'?^ A™"^ . . . A"^""]) - J) = 0. (55) 



A necessary stability condition for fractional order systems (21 ) to remain chaotic is keeping 
at least one eigenvalue A in the unstable region [62]. The number of saddle points and eigenval- 
ues for one-scroll, double-scroll and multi-scroll attractors was exactly described in work [63]. 
Assume that 3D chaotic system has only three equilibria. Therefore, if system has double-scroll 
attractor, it has two saddle points surrounded by scrolls and one additional saddle point. Sup- 
pose that the unstable eigenvalues of scroll saddle points are: Ai,2 = 0:1,2 ± j75i,2- The necessary 



condition to exhibit double-scroll attractor of system (21) is the eigenvalues A12 remaining in 



the unstable region [63]. The condition for commensurate derivatives order is 

g > -atan ( ^ ) , i = 1,2. (56) 



This condition can be used to determine the minimum order for which a nonlinear system can 
generate chaos [62]. 

Example 7. Let us investigate the Chen system with a double scroll attractor. The fractional 
order form of such system can be described as [66] 

D^-^Xi{t) = 35[x2(t) -xi(t)] 

Dl-°X2{t) = -7xi(t) -a;i(t)x3(t) + 28x2(t) 

D°-^X3(t) = xi(t)x2(t) -3x3(t) (57) 
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The system has three equihbrium at (0,0,0), (7.94,7.94,21) 
cobian matrix of the system evaluated at {xl,X2,x^) is: 

-35 35 



-7 



28 ~xl 
-3 



x; 



and (-7.94,-7.94,21). The Ja- 



(58) 



The two last equilibrium points are saddle points and surrounded by a chaotic double scroll 



attractor. For these two points, equation (55) becomes as follows: 

A^^ + 35X^^ + 3A^^ - 28A^^ + 105A^° - 21A^ + 4410 = 



(59) 



The characteristic equation (59) has unstable roots Ai^2 = 1.2928±0.2032j, |arg(Ai^2)| = 0.1560 



and therefore system (57) satisfy the necessary condition for exhibiting a double scroll attractor. 



Numerical simulation of the system (57) for initial conditions (—9, —5, 14) is depicted in Fig. ^ 




-20 30 



Xo(t) 



Figure 9: Double scroll attractor of Chen system (57) projected into 3D state space for 30 sec 



5 Conclusions 

In this paper we have presented the definitions for internal and external stability condition of 
certain class of the linear and nonlinear fractional order system of finite dimension given in 
state space, FODE or transfer function representation (polynomial). It is important to note 
that stability and asymptotic behavior of fractional order system is not exponential type [11] 
but it is in form of power law (a G i?), so called long memory behavior [36]. 

The results presented in this article are also applicable in robust stability investigation 
[24,47-49], stability of delayed system [15,22] and stability of discrete fractional order system 
[20,37]. Investigation of the fractional incommensurate order systems in state space, where 
space is deformed by various order of derivatives in various directions is still open. 
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